Lanczos Pseudospectral Propagation Method 
for Initial- Value Problems in Electrodynamics 

of Passive Media 



Andrei G. Borisov a ' x and Sergei V. Shabanov & ' 2 

a Laboratoire des Collisions Atomique et Moleculaires, UMR CNRS - Universite Paris-Sud 

8625, 91405 Orsay Cedex, France 

b Department of Mathematics, University of Florida, Gainesville, FL 32611, USA 

Abstract 

Maxwell's equations for electrodynamics of dispersive and absorptive (passive) me- 
dia are written in the form of the Schrodinger equation with a non-Hermitian Hamil- 
tonian. The Lanczos time-propagation scheme is modified to include non-Hermitian 
Hamiltonians and used, in combination with the Fourier pseudospectral method, to 
solve the initial-value problem. The time-domain algorithm developed is shown to be 
unconditionally stable. Variable time steps and/or variable computational costs per 
time step with error control are possible. The algorithm is applied to study transmis- 
sion and reflection properties of ionic crystal gratings with cylindric geometry in the 
infra-red range. 
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1. Introduction. There is a demand for fast time-domain solvers of the Maxwell's equa- 
tions to model dynamics of broad band electromagnetic pulses in dispersive and absorptive 
media. This is driven by various applications that include photonic devices, communications, 
and radar technology, amongst many others. Efficiency, accuracy, and stability are the key 
criteria of choosing a concrete algorithm for specific applications. In the past decades, pseu- 
dospectral methods of solving the initial value problem for differential equations have been 
under intensive study [1]. Because of their high efficiency and accuracy, they have replaced 
finite differencing approaches in many traditional applications as well as scientific simula- 
tions, e.g., in quantum chemistry [2]. Unconditionally stable pseudospectral algorithms are 
particularly attractive for numerical simulations. 

In the present paper we develop an unconditionally stable time-domain algorithm for solv- 
ing the initial value problem for Maxwell's equations in dispersive and absorptive (passive) 
media with sharp interfaces (discontinuities of medium parameters). It is a time-stepping 
algorithm that is based on the Hamiltonian formalism for electrodynamics of passive contin- 
uous media, the Lanczos propagation scheme [3, 4], and the Fourier pseudospectral method 
[5]. Apart from the unconditional stability, the algorithm has a dynamical control of accu- 
racy, which allows one to automatically optimize computational costs with error control at 
each time step. 

We apply the algorithm to the scattering of broad band electromagnetic pulses on grat- 
ings, the photonic devices that currently attract lots of attention because of their transmis- 
sion and reflection properties [6]. As for the passive medium, we choose an ionic crystal 
material. From the numerical point of view, the model of the dielectric permeability of 
such a material is rather representative and used in the vast number of applications. From 
the physical point of view, the interest to gratings and photonic crystals made of this kind 
of material is due two types of effects in interaction with electromagnetic radiation: The 
structural and polaritonic ones [7, 8]. We show that in the infrared range the reflection and 
transmission properties of ionic crystal gratings change significantly in narrow frequency 
ranges due to structural and polaritonic resonances. Structural resonances are associated 
with the existence of trapped (quasistationary) electromagnetic modes supported by the 
grating geometry (guided wave resonances) [9]. Polaritonic resonances are associated with 
dispersive properties of the material. Such resonances appear when the incident radiation 
can cause polaritonic excitations in the medium. From the macroscopic point of view, this 
occurs in the anomalous dispersion region of the dielectric constant. 

2. Basic equations. Maxwell's equations in passive media can be written in the form of 
the Schrodinger equation in which the wave function is a multidimensional column, composed 
of electromagnetic field components and the medium polarization, and the Hamiltonian is, 
in general, non-Hermitian when attenuation is present. The initial- value problem (the time 
evolution of an electromagnetic pulse) is then solved by finding the fundamental solution 
(the evolution operator kernel) for the Schrodinger equation. Here this idea is applied to 
the ionic crystal material whose dielectric properties at the frequency u are described by the 
dielectric constant 

e(u) = e QO + — ; — , 1 
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where e^o are constants, ujt is the resonant frequency, and r] is the attenuation. In particular, 
transmission and reflection properties of the periodic grating structure of circular parallel 
cylinders made of such a material are studied. 

Let P be a dispersive part of the total polarization vector of the medium. Then D = 
SooF, + P, where D and E are the electric induction and field, respectively. By using the 
Fourier transform, it is straightforward to deduce that P satisfies the second-order differential 
equation 



P + T]P + w£P 



(2) 



where the overdot denotes the time derivative, up 1 = (e — e^cu^/e^ if — £ oo is positive, 
otherwise, up 1 — > — in (2). Equation (2) must be solved with the zero initial conditions, 
p = p = o at t = 0. 

Define a set of auxiliary fields Qi j2 by P = y/e^u> p Qi / cut and Qi = u>tQ2- For non- 
magnetic media (/i — 1), the Maxwell's equations and (2) can be written as the Schrodinger 
equation, iip{t) = Hip{t), in which the wave function and the Hamiltonian are defined by 
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where c is the speed of light in the vacuum, and B is the magnetic field. A solution to the 
initial-value problem is given by ip(t) = exp(—itH)ip(0). Boundary conditions at medium 
interfaces are enforced dynamically, that is, parameters of the Hamiltonian H are allowed 
to be discontinuous functions of position. In particular, are set to one in the vacuum, 
and to some specific values in the medium in question (see Section 4). 



The norm of the wave function, \\i/j\\ 2 = J drip^ip, is proportional to the total electromag- 
netic energy of the wave packet [11, 12]. When the attenuation is not present, r] — 0, the 
Hamiltonian is Hermitian, W = H, relative to the conventional scalar product in the space 
of square integrable functions, and the norm (energy) is conserved. 

3. The algorithm. The Lanczos propagation scheme is based on a polynomial approx- 
imation of the short-time fundamental solution for the Schrodinger equation, ip(t + At) = 
exp(—iAtH)ip(t) [4, 10]. If the exponential is expanded into the Taylor series and the latter 
is truncated at the order 0(At n ), then the approximation of ip(t + At) belongs to the Krylov 
space K n = Span {ip(t), Hip(t), H n ~ l %p{t)}. Let P n be a projection operator of the original 
Hilbert space onto the Krylov space, P% = P n and = P n . Let P n ip = ip( n )- The exact 
solution of the initial-value problem is approximated by a solution of the corresponding ini- 
tial value problem in K n , ip {n) (t + At) = exp(-iAtH^)^ (n) (t), where # (n) = P n HP n . The 
accuracy of the approximation is of order 0(At n ). Hence, for practical needs, it is sufficient 
to choose n not so large (typically, n < 9). Then is a small matrix whose exponential is 
computed by direct diagonalization. If the Hamiltonian is Hermitian, then the corresponding 
matrix is symmetric tridiagonal in the Lanczos basis for K„ [3, 4]. Therefore the propaga- 
tion scheme is unitary, that is, the energy (norm) of the initial wave packet is preserved, 
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which also implies the unconditional stability of the algorithm. Another important feature 
of the Lanczos propagation scheme for Hermitian Hamiltonians is the dynamical control of 
accuracy [4, 10]. 

For non- Hermitian Hamiltonians, one possibility is to use the split method in combination 
with the Lanczos propagation scheme for the Hermitian part of the Hamiltonian [13] . While 
the unconditional stability is maintained, the accuracy is limited by the accuracy of the 
split approximation of the infinitesimal fundamental solution, typically, by 0(At 3 ). An 
alternative, mathematically sound, and more accurate procedure is based on the use of the 
dually orthonormal Lanczos basis in which retains the tridiagonal complex symmetric 
structure [3]. In the particular case of a complex symmetric H, the use of the dual Krylov 
space can be avoided. There is an orthogonal Lanczos basis for K n in which is also 
complex symmetric tridiagonal, but the orthogonality is now understood with respect to a 
new scalar product (without complex conjugation of vectors) [14]. However, in both the 
cases, the projector P n is no longer Hermitian. This has an unpleasant consequence. One 
can show that, while the accuracy remains of the order 0(At n ), the unconditional stability 
is typically lost. The algorithm is only conditionally stable. A more detailed study of such 
schemes will be given elsewhere. Here we focus on developing an unconditionally stable 
algorithm. For this purpose we construct an orthonormal basis for K„ itself, ignoring the 
dual Krylov space and making P n Hermitian. In this basis has a Hessenberg form, 

that is, it is upper-triangular with one extra non- zero lower superdiagonal [14]. A direct 
diagonalization of such a matrix is still not expensive for small n. The dynamical control of 
accuracy is also preserved. 

Let ipj, j = 0,1,. — 1, form an orthonormal basis for K n that is constructed as 
follows. Set 0o = i>(t) and ipo = 4>o/\\<j)o\\- Compute = (ipo,Hif) ) where (• , •) denotes 
the conventional scalar product in the Hilbert space of square integrable functions. For 
j = 1, 2, n — 1, compute 

h = ( 4 ) 

k=0 

h = h/UA , fc£-i = (ifc,ffifc-i) , (5) 

and, for k — 0, 1, j, compute 

h^ = (^,H^). (6) 

By construction, (ipj,ipk) = S jk . Let ip( n )(t + At) = P n ip(t + At) = Ylk=o c k(At)ip k where 
Cfc(At) = (ip k ,t/j(t + At)). Note that the basis functions ipj depend on t and so do Cj. For 
brevity of notations, the dependence on t of all quantities involving ipj is not explicitly shown 
in what follows. The expansion coefficients c k are united into a complex column c with n 
entries. It is straightforward to infer that c(At) satisfies the Schrodinger equation it = h^c 
with the initial condition c k (0) = S k o- The matrix hf^> has a Hessenberg form and is the 
projected Hamiltonian H^> in the basis constructed. Thus c(At) = exp(— iAth^)c(0). The 
exponential of the Hessenberg matrix is computed by direct diagonalization. 

To compute the basis functions, multiple actions of the Hamiltonian H on the current 



4 



wave function ip(t) are required. This is done by the Fourier pseudospectral method on a 
finite grid [5]. To suppress reflections of the wave packet from grid boundaries, absorbing 
boundary conditions are enforced by a conducting layer placed at the grid edges (see, e.g., 
[15]). The position dependence of the conductivity a is adjusted to suppress reflections 
with desired accuracy in the frequency range of interest. The Hamiltonian H is modified 
accordingly. In the upper left corner of H in (3), the function —inia is inserted instead of 
zero. 

Let us discuss the stability of the algorithm. In a time-stepping algorithm, an amplifica- 
tion matrix G(At) is defined by ip(t + At) = G(At)ip(t). The algorithm is unconditionally 
stable if HG^At)!! < const uniformly for all integers N > 0, At > and all other parameters 
characterizing the system [16]. The norm of an operator is defined by ||G|| = sup ||G'0||/||'0||- 
For any H the following decomposition holds: H = H + iV where H and V are Hermitian. 
Clearly, V is responsible for attenuation in any physically reasonable model of a passive 
medium. The norm of the initial- value problem solution decreases with time if V is negative 
semidefinite, that is, for any ip, (ip, Vip) < 0. Now we prove that the Lanczos algorithm with 
the Hessenberg projected Hamiltonian (6) is unconditionally stable, provided V of the total 
Hamiltonian is negative semidefinite. 

Observe that || exp(-iAtH) || < 1 for At > and all parameters of H for which V remains 
negative semidefinite. The amplification matrix has the form G(At) = exp(-iAtH^). 
Thanks to the Hermiticity of P„, it is sufficient to show that = P n VP n is negative 
semidefinite because the latter implies that HG^At)!! < HG^At)]^ < 1 uniformly for all 
integers N > 0, At > 0, and all parameters of H . For any if), the following chain of equalities 
holds, (ip,V^ip) = (ip , P n V Puip) = (?P(n), V^in)) < 0. In the first equality, the definition 
of has been used, in the second one, the Hermiticity of the projection operator has 
been invoked, and the final inequality is valid since V is negative semidefinite. The proof is 
completed. 

The most expensive operation of the algorithm is the action of the Hamiltonian on the 
wave function. Bearing also in mind the accumulation of round-off errors, when computing 
powers of H with broad spectrum (cf. [3]), this implies that the dimension of the Krylov 
space has to be as small as possible at each time step while controlling the approximation 
error. In the case of Hermitian Hamiltonians, the n needed at a specified time step can be 
deduced from the condition [4] that |c n _i(At)| 2 < e where e is a small number. This condition 
is based on the fact that c n _i(At) ~ 0(At n_1 ). Note that c n _i determines the weight of 
H n - 1 ip(t) in the Taylor expansion of ip{t + At). In our algorithm for a non- Hermitian H, 
the time evolution of the vector c is generated by a Hessenberg matrix h^ n \ By examining 
the Taylor expansion of the exponential in exp(— iAth^)c(0) it is easy to convince oneself 
that c„_i(At) ~ 0(At n_1 ) remains valid thanks to Cj(0) = 5j . Thus, the same dynamical 
accuracy control can be used. In particular, in our simulations the dimension of the Krylov 
space at each time step is determined by |c n _ 3 | 2 + |c„_ 2 | 2 + |c n _i | 2 < e at a fixed At to control 
weights of three highest Krylov vectors in the approximate solution, and e ~ 10~ 14 with At 
being fixed. Note, however, that both the propagation parameters n and At can be varied 
at each time step to minimize computational costs at the given accuracy level. 
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4. Ionic crystal gratings. The algorithm has been applied to simulate the scattering 
of broad band electromagnetic (laser) pulses on a grating structure consisting of circular 
parallel ionic crystal cylinders periodically arranged in vacuum. Our primary interest is to 
study the effect of trapped modes (guided wave resonances) and polaritonic excitations on 
transmission and reflection properties of the grating in the infrared range. The dielectric 
function of the ionic crystal material is approximated by the single oscillator model (1). 
Following the work [7], we choose the parameters representative for the beryllium oxide: 
£oo = 2.99, e = 6.6, ojt = 87.0 meV, and the damping rj = 11.51 meV. The packing density 
R/D g = 0.1, where R is the radius of cylinders and D g is the grating period, has been kept 
fixed in simulations. The cylinders are set parallel to the y axis. The structure is periodic 
along the x axis, and the z direction is transverse to the grating. A Gaussian wave packet 
propagating along the z axis is used as an initial configuration. It is linearly polarized with 
the electric field oriented along the y axis, i.e., parallel to the cylinders (the so called TM 
polarization). The frequency resolved transmission and reflection coefficients are obtained 
via the time-to-frequency Fourier transform of the signal on "virtual detectors" placed at 
some distance in front and behind the periodic structure [17]. The zero diffraction mode 
is studied here for wavelengths A > D g so that reflected and transmitted beams propagate 
along the z-axis. Similar to our previous works [11, 13] we use a change of variables in both 
x (x = an d z (z — /2(C)) coordinates to enhance the sampling efficiency in the vicinity 

of medium interfaces. A typical size of the mesh corresponds to —17D g < z < 15D g , and 
—0.5D g < x < 0.5D g with, respectively, 384 and 64 knots. Note that, because of the variable 
change, a uniform mesh in the auxiliary coordinates (£, () corresponds to a non-uniform mesh 
in the physical (x, z) space. 

Two types of resonances are expected for the gratings studied here. Structure resonances 
are associated with the existence of guided wave modes [9, 13]. They are characteristic for 
periodic dielectric gratings and, in the absence of losses, lead to the 100% reflection within 
a narrow frequency interval(s) for wavelengths A ~ D g . The second type of resonances arise 
because of polariton excitations for wavelengths A ~ Dt = 2-kc/uot- Calculations have been 
done for different values of D g so that the polaritonic excitation can be tuned through out 
the wavelength range of interest (\/D g > 1) by changing the ratio Dt/ D g . 

In Fig. 1 we show the results obtained for the transmission (blue curves) and reflection 
(red curves) coefficients for the beryllium oxide gratings characterized by the period D g such 
that Dt/D 9 = 0.5, 2.5, and 4 as indicated in the figure. The results are presented as a 
function of the radiation wavelength measured in units of the grating period. Note that the 
logarithmic scale is used for the horizontal axis in order to improve the resolution at small 
wavelengths. Consider first the following two limiting cases. According to (1), for short 
wavelengths A <C Dt {D g <C Dt), the medium behaves dielectric with £ ^ £ m . In 

the long wavelength limit A ^> Dt {D g ^> Dt), the medium responds as a dielectric material 
characterized by e ~ sq. In Fig. 1 the dashed and solid black curves represent the reflection 
coefficient of the grating made of a lossless, non-dispersive dielectric with e = and e = Eo, 
respectively. In agreement with the previously published results [13], the reflection coefficient 
in these two cases reaches 1 within a narrow frequency range for A ~ D g . This resonant 
pattern is associated with the so-called Wood anomalies [18], and can be explained by the 



6 



existence of trapped modes or guided wave resonances [9, 11]. The width of the resonances 
is determined by the lifetime of the corresponding quasi-stationary trapped mode which is 
a standing wave along the x axis and is excited by the incoming wave. The width increases 
with e while the resonant wavelength gets redshifted, which explains the difference between 
the dashed and solid black curves (e > £qo) . 

Now we turn to the discussion of the effects due to dispersive properties of the ionic crystal 
material. For D T /D g = 0.5, the resonant excitation of polaritons is impossible within the 
range of wavelengths of interest, and the dielectric constant is close to Eq. The result for 
the reflection coefficient in this case is similar to the data shown by the black solid curve. 
However, there is an essential difference as compared to the case of a lossless, non- dispersive 
dielectric grating. Indeed, for a lossless medium the sum of the reflection and transmission 
coefficients must be one, which is not the case for the beryllium oxide model because of the 
damping (the dashed-dotted green curve). The maximal loss of energy corresponds to the 
resonant wavelength. It is easily understood because the trapped mode remains in contact 
with the material much longer than the main pulse, and, therefore, can dissipate more energy. 

For Dx/Dg = 2.5, two resonances emerge leading to the enhanced reflection within the 
corresponding frequency ranges. The one at X/D g ~ 2.5, i.e., A ~ D?, is associated with 
polaritonic excitations of the ionic crystal. The resonance at A ~ D g is a structure resonance. 
As follows from (1), the dielectric constant in this case approaches e^ for small wavelengths 
A ~ D g . Then the width and position of the structure resonance are close to the data given 
by the dashed black curve. The imaginary part of the dielectric function is large enough 
through out the entire wavelength range to produce a substantial energy loss at both the 
resonances. 

Finally, for Dt/D 9 = 4 the polariton excitation appears at A ~ 4D g and the two reso- 
nances are well separated. The structure resonance at A ~ D g closely matches the result for 
a lossless, non-dispersive dielectric grating characterized by e = e^. Observe that the reflec- 
tion coefficient is close to 1 in this case and the energy loss is small because the imaginary 
part of e{uj) is small far from uj = ojt- 

Figs 2 and 3 show the reflection and transmission coefficients of the grating as functions 
of the incident radiation wavelength and the grating period D g . The polariton resonance 
wavelength D T = 2itc/uj t and the packing density R/D g are kept fixed. The results for the 
two limiting cases in which e = Eq and e = £oo are represented by the left most and right 
most colored columns, respectively. The resonance pattern of the system is clearly visible, in 
particular, the transformation of the structure resonance at e = £q into the polaritonic one. 
Thus, by increasing the ratio Dt/ D g , the "broad" structure resonance associated with e = £o 
is turned into the polaritonic resonance and follows the diagonal of the plot (X/D g = D T /D g ). 
At the same time, starting approximately with D T /D g = 2, the "narrow" structure resonance 
associated with £ = £00 emerges and fully develops for Dt/ D g = 4. 

Finally we would like to show the sensitivity of the results to the attenuation in the 
present system. In Fig. 4, the transmission and reflection coefficients are presented for two 
different choices of the attenuation rj in (1). The geometry of the grating structure is set 
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by D T /D g = 2.5. The upper panel of the figure corresponds to 77 = 11.51 meV as used 
throughout this paper. The lower panel of the figure corresponds to the damping reduced 
by the factor of 20 : rj — > rj/20. Overall features are qualitatively the same in both the cases. 
Thus, both the structure resonance at A ~ D g and polaritonic resonance at A ~ 2.5D g 
are present, accompanied by the reduced transmission and enhanced reflection. As the 
wavelength increases from the structure resonance, the reflectivity of the grating drops to 
zero. Its subsequent onset for A > 1.683D 9 is linked with the metal-type behavior of the ionic 
crystal (e becomes negative). The characteristic frequency for the "metallization" can be 
deduced from the Lyddane-Sachs- Teller relation u L = ut^Eq/ 'e^, leading to A^ = 1.683D g 
for Ay = 2.5D g . Despite these common features, the reduction of the attenuation leads to 
essential changes. In contrast to the upper panel of the figure, for 77 — > 77 / 20 the transmission 
coefficient reaches nearly at both the resonances, and the reflection coefficient is close to 1. 
Moreover, new structures appear in the polaritonic resonance for A = Dt, i.e., as e changes 
from large negative to large positive values. These structures are completely washed out for 
the medium with large damping. This result indicates the importance of accurate modelling 
of losses in polaritonic media in order to make reliable predictions of transmission and 
reflection properties of grating structures and photonic crystals. 

5. Conclusions. We have developed an unconditionally stable (time-domain) algorithm 
for initial value problems in electrodynamics of inhomogeneous, dispersive, and absorptive 
media. The method is based on the three essential ingredients: (i) the Hamiltonian formal- 
ism in electrodynamics of passive media, (ii) the Lanczos propagation scheme, modified to 
account for attenuation, and (ra) the Fourier pseudospectral method on non-uniform grids 
induced by change of variables to enhance the sampling efficiency in the vicinity of sharp 
inhomogeneities of the medium. Apart from the unconditional stability, the algorithm al- 
lows for a dynamical accuracy control, meaning that the two propagation parameters, the 
dimension of the Krylov space and the time step, may automatically be adjusted to minimize 
computational costs in due course of simulations, while controlling error. 

The algorithm has been tested by simulating the scattering of infrared electromagnetic 
pulses on periodic gratings composed of parallel cylinders that are made of the ionic crystal 
material. The Lorentz model describing dielectric properties of such a material is rather 
representative and used to model a vast variety of dielectric materials. Our results demon- 
strate the role of structure (or guided wave) resonances and polaritonic excitations for the 
transmission and reflection properties of grating structures. The results are also shown to 
be sensitive to the attenuation of polaritonic media. 
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Figure captions 

Fig. 1. Calculated zero-order reflection (red curves) and transmission (blue curves) 
coefficients for the ionic crystal grating described in the text. The results are presented as 
a function of the incident radiation wavelength measured in units of the grating period D g . 
Different panels of the figure correspond to different values of the grating period as compared 
to the resonance wavelength for the polaritonic excitation of the material, Dt = 2nc/uT- 
The dashed and solid black curves represent the reflection coefficient calculated for the 
grating made of a lossless, non- dispersive dielectric characterized by e = and e = e , 
respectively. The sum of the reflection and transmission coefficients is shown as the dashed- 
dotted green curve. Its deviation from 1 represents the electromagnetic energy loss because 
of the attenuation. 

Fig. 2. The zero-order transmission coefficient for the ionic crystal grating described 
in the text as a function of the incident radiation wavelength and the grating period. The 
horizontal axis represents the ratio Dt/D 9 of the resonance wavelength for the polaritonic 
excitation of the material Dt = 2-kc/ut and the grating period D g . The vertical axis 
represents the incident radiation wavelength A measured in units of D g . Color codes used 
for the plot are shown in the inset. 

Fig. 3. The zero-order reflection coefficient for the ionic crystal grating described in 
the text as a function of the incident radiation wavelength and the grating period. The 
horizontal axis represents the ratio D T / D g of the resonance wavelength for the polaritonic 
excitation of the material Dt = 2ttc/ujt and the grating period D g . The vertical axis 
represents the incident radiation wavelength A measured in units of D g . Color codes used 
for the plot are shown in the inset. 

Fig. 4. Calculated zero-order reflection (red curves) and transmission (dashed blue 
curves) coefficients for the ionic crystal grating. The sum of the reflection and transmission 
coefficients is shown as the dashed-dotted green curve. The geometry of the grating structure 
is set by D T /D g = 2.5. The upper panel of the figure corresponds to the attenuation 
rj = 11.51 meV as used throughout the paper. The lower panel of the figure corresponds 
to the damping reduced by the factor of 20 : rj — > r]/20. The vertical black line defines the 
resonant wavelength A = D T = 2-nc/u T . 
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